Exhausted CD8 T cells are associated with ICB responsiveness in breast cancer

Background

The single-cell dataset generated by Bassez et al. (2021) Nature Medicine contains single-cell data from breast carcinoma biopsies from 29 patients taken immediately before (pre-treatment) and 9 +/- 2 days after (on-treatment) receiving anti-PD1 immunotherapy. Lacking an objective measure of clinical benefit following therapy, in this study the authors use T cell expansion following treatment as a proxy for response.

Full Data Original Source

Case Study Dataset Internal link

Goal

To compare the T cell subtype composition between responders and non-responders using ProjecTILs and a reference TIL atlas

R Environment

# renv::activate()
renv::restore()
# remotes::install_github('carmonalab/ProjecTILs', ref='1.0.0')
# remotes::install_version('Seurat', version = '4.0.1')
library(ProjecTILs)
library(patchwork)
library(ggplot2)
library(reshape2)
library(parallelly)
options(parallelly.makeNodePSOCK.setup_strategy = "sequential")  # workaround to issue https://github.com/HenrikBengtsson/future/issues/511

scRNA-seq data preparation

ddir <- "input/Bassez_breast_data"
if (!file.exists(ddir)) {
    dir.create(ddir)
    dataUrl <- "https://drive.switch.ch/index.php/s/7jdqnvPZuJrriZB/download"
    download.file(dataUrl, paste0(ddir, "/tmp.zip"))
    unzip(paste0(ddir, "/tmp.zip"), exdir = ddir)
    file.remove(paste0(ddir, "/tmp.zip"))
}



# Count matrices
f1 <- sprintf("%s/1864-counts_tcell_cohort1.rds", ddir)

cohort1 <- readRDS(f1)
dim(cohort1)
[1] 25288 53382
# Metadata
meta1 <- read.csv(sprintf("%s/1870-BIOKEY_metaData_tcells_cohort1_web.csv", ddir))
rownames(meta1) <- meta1$Cell

head(meta1)
                                                             Cell nCount_RNA
BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1 BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1       1252
BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1 BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1       1869
BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1 BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1       1000
BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1 BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1       1288
BIOKEY_13_Pre_AAAGATGGTTTGACAC-1 BIOKEY_13_Pre_AAAGATGGTTTGACAC-1       2056
BIOKEY_13_Pre_AAAGATGTCAACGCTA-1 BIOKEY_13_Pre_AAAGATGTCAACGCTA-1       1224
                                 nFeature_RNA patient_id timepoint expansion
BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1          700  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1          995  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1          627  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1          681  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAAGATGGTTTGACAC-1          789  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAAGATGTCAACGCTA-1          634  BIOKEY_13       Pre       n/a
                                 BC_type cellSubType          cohort
BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1   HER2+         gdT treatment_naive
BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1   HER2+     NK_REST treatment_naive
BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1   HER2+       CD4_N treatment_naive
BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1   HER2+      CD8_EM treatment_naive
BIOKEY_13_Pre_AAAGATGGTTTGACAC-1   HER2+       CD8_N treatment_naive
BIOKEY_13_Pre_AAAGATGTCAACGCTA-1   HER2+      CD8_RM treatment_naive
data.seurat <- CreateSeuratObject(cohort1, project = "Cohort1_IT", meta.data = meta1)

Subset pre-treatment biopsies

table(data.seurat$timepoint)

   On   Pre 
27528 25854 
data.seurat <- subset(data.seurat, subset = timepoint == "Pre")

Subset T cells according to annotation by the authors

table(data.seurat$cellSubType)

               CD4_EM                CD4_EX  CD4_EX_Proliferating 
                 4810                  1707                   107 
                CD4_N               CD4_REG CD4_REG_Proliferating 
                 3007                  2861                    69 
               CD8_EM              CD8_EMRA                CD8_EX 
                 5443                   290                  2286 
 CD8_EX_Proliferating                 CD8_N                CD8_RM 
                  340                   354                  1963 
                  gdT               NK_CYTO               NK_REST 
                 1000                   230                   938 
           Vg9Vd2_gdT 
                  449 
data.seurat <- subset(data.seurat, subset = cellSubType %in% c("NK_REST", "Vg9Vd2_gdT",
    "gdT", "NK_CYTO"), invert = T)

Basic QC, remove outliers in # counts, # detected genes, and mitochondrial and ribosomal content

percent.ribo.dv <- PercentageFeatureSet(data.seurat, pattern = "^RP[LS]")
percent.mito.dv <- PercentageFeatureSet(data.seurat, pattern = "^MT-")

data.seurat <- AddMetaData(data.seurat, metadata = percent.ribo.dv, col.name = "percent.ribo")
data.seurat <- AddMetaData(data.seurat, metadata = percent.mito.dv, col.name = "percent.mito")

Idents(data.seurat) <- "patient_id"
VlnPlot(data.seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.ribo", "percent.mito"),
    ncol = 2, pt.size = 0)

data.seurat <- subset(data.seurat, subset = nFeature_RNA > 500 & nFeature_RNA < 6000 &
    nCount_RNA > 600 & nCount_RNA < 25000 & percent.ribo < 50 & percent.mito < 15)

Remove samples that are too small and downsample large ones, in order to have similar contribution from all patients/samples contribute. Downsampling large samples also speeds up downstream computations.

minCells <- 500
ds <- 1000

Idents(data.seurat) <- "patient_id"
table(data.seurat$patient_id)

 BIOKEY_1 BIOKEY_10 BIOKEY_11 BIOKEY_12 BIOKEY_13 BIOKEY_14 BIOKEY_15 BIOKEY_16 
     3014      1006       543      2511      1220       920      1026      2070 
BIOKEY_17 BIOKEY_18 BIOKEY_19  BIOKEY_2 BIOKEY_20 BIOKEY_21 BIOKEY_22 BIOKEY_23 
       96       182      1430       945       109       194        68        35 
BIOKEY_24 BIOKEY_25 BIOKEY_26 BIOKEY_27 BIOKEY_28 BIOKEY_29  BIOKEY_3 BIOKEY_30 
      339        95        85       594       575       162       264        61 
BIOKEY_31  BIOKEY_4  BIOKEY_5  BIOKEY_6  BIOKEY_7  BIOKEY_8  BIOKEY_9 
      346      2282      1688       606        63       135       168 
largeSamples <- names(table(data.seurat$patient_id)[table(data.seurat$patient_id) >=
    minCells])
data.seurat <- subset(data.seurat, idents = largeSamples)
table(data.seurat$patient_id)

 BIOKEY_1 BIOKEY_10 BIOKEY_11 BIOKEY_12 BIOKEY_13 BIOKEY_14 BIOKEY_15 BIOKEY_16 
     3014      1006       543      2511      1220       920      1026      2070 
BIOKEY_19  BIOKEY_2 BIOKEY_27 BIOKEY_28  BIOKEY_4  BIOKEY_5  BIOKEY_6 
     1430       945       594       575      2282      1688       606 
data.seurat <- subset(data.seurat, cells = WhichCells(data.seurat, downsample = ds))
table(data.seurat$patient_id)

 BIOKEY_1 BIOKEY_10 BIOKEY_11 BIOKEY_12 BIOKEY_13 BIOKEY_14 BIOKEY_15 BIOKEY_16 
     1000      1000       543      1000      1000       920      1000      1000 
BIOKEY_19  BIOKEY_2 BIOKEY_27 BIOKEY_28  BIOKEY_4  BIOKEY_5  BIOKEY_6 
     1000       945       594       575      1000      1000       606 

ProjecTILs analysis

Here we will project the query data onto a reference TIL atlas to interpret cell states in the context of reference T cell subtypes.

For optimal batch-effect correction, we prefer projecting each patient/batch separately (although results tend to be very similar if projecting all together)

Here we use filter.cells=F to save processing time, because we rely on the authors’ cell type annotation to subset T cells.

For large samples make.projection can take several minutes. Set the number of parallel jobs with ncores according to your computational resources (you can check future::availableCores(), although parallelization is usually limited by available RAM)

ncores = 5
ref <- load.reference.map()  # by default, murine reference TIL atlas v1.0
[1] "Loading Default Reference Atlas..."
[1] "/Users/admin/gitHub/ProjecTILs_CaseStudies/ref_TILAtlas_mouse_v1.rds"
[1] "Loaded Reference map ref_TILAtlas_mouse_v1"
data.split <- SplitObject(data.seurat, split.by = "patient_id")
query.projected <- make.projection(data.split, ref, filter.cells = F, ncores = ncores)
[1] "Using assay RNA for BIOKEY_13"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_13 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_10"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_10 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_16"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_16 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_14"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_14 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_19"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_19 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_28"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_28 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_15"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_15 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_5"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_5 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_12"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_12 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_1"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_1 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_4"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_4 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_11"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_11 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_2"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_2 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_6"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_6 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"
[1] "Using assay RNA for BIOKEY_27"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning BIOKEY_27 to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"

See query cells for two selected the samples, projected onto the reference atlas

plot.projection(ref, query = query.projected$BIOKEY_13) + ggtitle("Projection sample 1")

plot.projection(ref, query = query.projected$BIOKEY_10) + ggtitle("Projection sample 2")

For simplified downstream analysis, merge back ProjecTILs results list of objects into a single Seurat object

query.projected <- lapply(query.projected, function(x) {
    cellstate.predict(ref = ref, query = x, reduction = "umap", ndim = 2)
})

query.projected.merged <- suppressMessages(Reduce(ProjecTILs:::merge.Seurat.embeddings,
    query.projected))

We visually inspect now how the gene profiles of projected cells match those of the reference, for each T cell subtype.

A query dataset might contain a subtype that is not well represented in the atlas. That would requiere splitting heterogeneous populations after projection by higher-dimensions (e.g. Find Discriminant Dimensions ) In this case, radar/spider profiles look fine, meaning that in general terms, current reference atlas is capturing the main features of the query data:

plot.states.radar(ref, query = query.projected.merged, min.cells = 30, genes4radar = c("Foxp3",
    "Cd4", "Cd8a", "Tcf7", "Ccr7", "Sell", "Il7r", "Gzmb", "Gzmk", "Pdcd1", "Lag3",
    "Havcr2", "Tox", "Entpd1", "Cxcr3", "Cxcr5", "Ifng", "Cd200", "Xcl1", "Itgae"))  #,'Cxcl13',

We can also verify that gene profiles are consistent across samples/patients, i.e. batch effects have been mitigated

plot.states.radar(ref, query = query.projected[1:5], min.cells = 30, genes4radar = c("Foxp3",
    "Cd4", "Cd8a", "Tcf7", "Ccr7", "Sell", "Il7r", "Gzmb", "Gzmk", "Pdcd1", "Lag3",
    "Havcr2", "Tox", "Entpd1", "Cxcr3", "Ifng", "Cd200", "Xcl1", "Itgae"))

Compare responders vs. non-responders (based on clonal expansion) - pre-therapy

# Split by expansion (defined by authors)
query.sub.byExp <- subset(query.projected.merged, subset = expansion %in% c("E",
    "NE"))
query.sub.byExp <- SplitObject(query.sub.byExp, split.by = "expansion")

n_e <- length(table(query.sub.byExp$E$patient_id))
n_ne <- length(table(query.sub.byExp$NE$patient_id))

a <- plot.projection(ref, query = query.sub.byExp$E, linesize = 0.5, pointsize = 0.2) +
    ggtitle(paste0("Responders; n=", n_e))
b <- plot.projection(ref, query = query.sub.byExp$NE, linesize = 0.5, pointsize = 0.2) +
    ggtitle(paste0("Non-Responders; n=", n_ne))

a

b

c <- plot.statepred.composition(ref, query.sub.byExp$E, metric = "percent") + ylim(0,
    30) + ggtitle(paste0("Responders; n=", n_e))
d <- plot.statepred.composition(ref, query.sub.byExp$NE, metric = "percent") + ylim(0,
    30) + ggtitle(paste0("Non-Responders; n=", n_ne))
c | d

The analysis shows that ICB-responding tumors have a higher proportion of exhausted (Tex) and precursor-exhausted (Tpex) CD8 TILs, as well as a lower ratio of Treg/Tex.

Calculate and visualize fold-changes of TIL composition in responders vs non-responders

which.types <- table(query.projected.merged$functional.cluster) > 20  # disregard very small populations
stateColors_func <- c("#edbe2a", "#A58AFF", "#53B400", "#F8766D", "#00B6EB", "#d1cfcc",
    "#FF0000", "#87f6a5", "#e812dd")
states_all <- levels(ref$functional.cluster)
names(stateColors_func) <- states_all
cols_use <- stateColors_func[names(which.types)][which.types]

query.projected.merged$functional.cluster <- factor(query.projected.merged$functional.cluster,
    levels = states_all)
query.list <- SplitObject(query.projected.merged, split.by = "expansion")

norm.c <- table(query.list[["NE"]]$functional.cluster)/sum(table(query.list[["NE"]]$functional.cluster))
norm.q <- table(query.list[["E"]]$functional.cluster)/sum(table(query.list[["E"]]$functional.cluster))

foldchange <- norm.q[which.types]/norm.c[which.types]
foldchange <- sort(foldchange, decreasing = T)

tb.m <- melt(foldchange)
colnames(tb.m) <- c("Cell_state", "Fold_change")
pll <- list()
ggplot(tb.m, aes(x = Cell_state, y = Fold_change, fill = Cell_state)) + geom_bar(stat = "identity") +
    scale_fill_manual(values = cols_use) + geom_hline(yintercept = 1) + scale_y_continuous(trans = "log2") +
    theme(axis.text.x = element_blank(), legend.position = "left") + ggtitle("Responder vs. Non-responder")

Tumors in ICB-responders have a ~4 fold-change increase in Tex and a ~2 fold-change increase in both Tpex and Tregs.

Discussion

Our automated ProjecTILs analysis of this patient cohort indicate that compared to non-responders, early responders to PD-1 blockade (in terms of T cell expansion) have higher relative amounts of CD8 Tex and Tregs infiltrating their tumors, and a lower Treg to T(p)ex ratio, in agreement with previous studies (Daud et al. 2016; Thommen et al. 2018; Kumagai et al. 2020; Bassez et al. 2020). This supports the notion that patients with a pre-existing tumor-specific CD8 T cell response in their tumors are more likely to respond to PD-1 blockade therapy.

Because tumoricidal CD8 Tex are derived from Tpex (Siddiqui et al. 2019; Miller et al. 2019), the presence of CD8 Tex seems to be a good proxy for the presence of smaller populations of quiescent Tpex, that have the ability to proliferate and differentiate into Tex cells following PD-1 blockade (Siddiqui et al. 2019; Miller et al. 2019). Indeed we observed a small population of Tpex-like cells, with increased levels R vs RN. Because there are very few Tpex their gene profile might be less robust than others, but compared to Tex they display higher levels of Tpex-associated genes, such as XCL1 and CD200, and lower levels of terminal-exhaustion markers such as HAVCR2 and ENTPD1.

In their study, Bassez et al. also observed increased levels in R vs NR of a population they referred to as exhausted CD4 T cells. Here we observed also a trend for increased levels of follicular-helper CD4 T cells, that display exhaustion features (e.g. TOX, PDCD1, ENTPD1 expression), although the most prominent enrichment in ICB-responders is for exhausted CD8 T cells.

Further reading

Dataset original publication - Bassez et al

ProjecTILs case studies - INDEX - Repository

The ProjecTILs method Andreatta et. al (2021) Nat. Comm. and code

References

  • Bassez, A., Vos, H., Van Dyck, L. et al. A single-cell map of intratumoral changes during anti-PD1 treatment of patients with breast cancer. Nat Med 27, 820–832 (2021). https://doi.org/10.1038/s41591-021-01323-8

  • Daud, A.I., Loo, K., Pauli, M.L., Sanchez-Rodriguez, R., Sandoval, P.M., Taravati, K., Tsai, K., Nosrati, A., Nardo, L., Alvarado, M.D., Algazi, A.P., Pampaloni, M.H., Lobach, I.V., Hwang, J., Pierce, R.H., Gratz, I.K., Krummel, M.F., Rosenblum, M.D., 2016. Tumor immune profiling predicts response to anti-PD-1 therapy in human melanoma. J. Clin. Invest. 126, 3447–3452. https://doi.org/10.1172/JCI87324

  • Gros, A., Robbins, P.F., Yao, X., Li, Y.F., Turcotte, S., Tran, E., Wunderlich, J.R., Mixon, A., Farid, S., Dudley, M.E., Hanada, K., Almeida, J.R., Darko, S., Douek, D.C., Yang, J.C., Rosenberg, S. a, 2014. PD-1 identifies the patient-specific in filtrating human tumors. J. Clin. Invest. 124, 2246–59. https://doi.org/10.1172/JCI73639.2246

  • Miller, B.C., Sen, D.R., Al Abosy, R., Bi, K., Virkud, Y.V., LaFleur, M.W., Yates, K.B., Lako, A., Felt, K., Naik, G.S., Manos, M., Gjini, E., Kuchroo, J.R., Ishizuka, J.J., Collier, J.L., Griffin, G.K., Maleri, S., Comstock, D.E., Weiss, S.A., Brown, F.D., Panda, A., Zimmer, M.D., Manguso, R.T., Hodi, F.S., Rodig, S.J., Sharpe, A.H., Haining, W.N., 2019. Subsets of exhausted CD8+ T cells differentially mediate tumor control and respond to checkpoint blockade. Nat. Immunol. 20, 326–336. https://doi.org/10.1038/s41590-019-0312-6

  • Siddiqui, I., Schaeuble, K., Chennupati, V., Fuertes Marraco, S.A., Calderon-Copete, S., Pais Ferreira, D., Carmona, S.J., Scarpellino, L., Gfeller, D., Pradervand, S., Luther, S.A., Speiser, D.E., Held, W., 2019. Intratumoral Tcf1+PD-1+CD8+ T Cells with Stem-like Properties Promote Tumor Control in Response to Vaccination and Checkpoint Blockade Immunotherapy. Immunity 50, 195–211.e10. https://doi.org/10.1016/J.IMMUNI.2018.12.021

  • Thommen, D.S., Koelzer, V.H., Herzig, P., Roller, A., Trefny, M., Dimeloe, S., Kiialainen, A., Hanhart, J., Schill, C., Hess, C., Prince, S.S., Wiese, M., Lardinois, D., Ho, P.C., Klein, C., Karanikas, V., Mertz, K.D., Schumacher, T.N., Zippelius, A., 2018. A transcriptionally and functionally distinct pd-1 + cd8 + t cell pool with predictive potential in non-small-cell lung cancer treated with pd-1 blockade. Nat. Med. 24, 994. https://doi.org/10.1038/s41591-018-0057-z

  • Kumagai, S., Togashi, Y., Kamada, T. et al. The PD-1 expression balance between effector and regulatory T cells predicts the clinical efficacy of PD-1 blockade therapies. Nat Immunol 21, 1346–1358 (2020). https://doi.org/10.1038/s41590-020-0769-3